#include "rayleigh_singlescat.h"
#include <math.h>

RayleighSingleScat::RayleighSingleScat(QString layerName, QObject *parent) :
    SingleScatLayer(layerName, parent)
{
}

void RayleighSingleScat::solarScatCoeff(int iThetaS, int iPhiS, linearPolar *lightOut)
{
    double phiS = iPhiS * M_PI / N_PI;
    double phiI = m_solarAzi;
    double thetaS = (iThetaS > 0) ?
		(iThetaS*M_PI/N_PI + ANGLE_STEP/2) :
		(M_PI + iThetaS*M_PI/N_PI -ANGLE_STEP/2);

    double thetaI = m_solarZen;

    double lightIncident[3];
    lightIncident[0] = m_solarIntensity ;
    lightIncident[1] = 0;
    lightIncident[2] = 0;

    lightOut->i = (3*lightIncident[U]*Sin(phiI - phiS)*Sin(thetaS)*(Cos(thetaS)*Sin(thetaI) - Cos(phiI - phiS)*Cos(thetaI)*Sin(thetaS)))/2. + (3*lightIncident[I]*(1 + Power(Cos(thetaI)*Cos(thetaS) + Cos(phiI - phiS)*Sin(thetaI)*Sin(thetaS),2)))/4. - (3*lightIncident[Q]*(2*Power(Cos(thetaS),2)*Power(Sin(thetaI),2) + 2*Power(Cos(phiI - phiS),2)*Power(Cos(thetaI),2)*Power(Sin(thetaS),2) - 2*Power(Sin(phiI - phiS),2)*Power(Sin(thetaS),2) - Cos(phiI - phiS)*Sin(2*thetaI)*Sin(2*thetaS)))/8.;

    lightOut->q = (3*lightIncident[U]*Sin(phiI - phiS)*(Cos(phiI - phiS)*Cos(thetaI)*(3 + Cos(2*thetaS)) + Sin(thetaI)*Sin(2*thetaS)))/4. - (3*lightIncident[I]*(2*Power(Cos(phiI - phiS),2)*Power(Cos(thetaS),2)*Power(Sin(thetaI),2) - 2*Power(Sin(phiI - phiS),2)*Power(Sin(thetaI),2) + 2*Power(Cos(thetaI),2)*Power(Sin(thetaS),2) - Cos(phiI - phiS)*Sin(2*thetaI)*Sin(2*thetaS)))/8. + (3*lightIncident[Q]*(Cos(2*phiI)*Cos(2*phiS)*(3 + Cos(2*thetaI))*(3 + Cos(2*thetaS)) + (3 + Cos(2*thetaI))*(3 + Cos(2*thetaS))*Sin(2*phiI)*Sin(2*phiS) + 12*Power(Sin(thetaI),2)*Power(Sin(thetaS),2) + 4*Cos(phiI - phiS)*Sin(2*thetaI)*Sin(2*thetaS)))/32.;

    lightOut->u = (3*lightIncident[I]*Sin(phiI - phiS)*Sin(thetaI)*(Cos(phiI - phiS)*Cos(thetaS)*Sin(thetaI) - Cos(thetaI)*Sin(thetaS)))/2. + (3*lightIncident[U]*(Cos(2*(phiI - phiS))*Cos(thetaI)*Cos(thetaS) + Cos(phiI - phiS)*Sin(thetaI)*Sin(thetaS)))/2. - (3*lightIncident[Q]*Sin(phiI - phiS)*(Cos(phiI - phiS)*(3 + Cos(2*thetaI))*Cos(thetaS) + Sin(2*thetaI)*Sin(thetaS)))/4.;

//    lightOut->i *= (1/(4*M_PI));
//    lightOut->q *= (1/(4*M_PI));
//    lightOut->u *= (1/(4*M_PI));
    lightOut->i /= 4;
    lightOut->q /= 4;
    lightOut->u /= 4;
}
